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ABSTRACT 

PSR B 1259-63/LS 2883 is a binary system in which a 48-ms pulsar orbits around 
a Be star in a high eccentric orbit with a long orbital period of about 3.4 yr. It is 
special for having asymmetric two-peak profiles in both the X-ray and the TeV light 
curves. Recently, an unexpected GeV flare was detected by Fermi gamma-ray obser- 
vatory several weeks after the last periastron passage. In this paper, we show that this 
observed GeV flare could be produced by the Doppler-boosted synchrotron emission 
in the bow shock tail. An anisotropic pulsar wind model, which mainly affects the 
energy flux injection to the termination shock in different orbital phase, is also used in 
this paper, and we find that the anisotropy in the pulsar wind can play a significant role 
in producing the asymmetric two-peak profiles in both X-ray and TeV light curves. 
The X-ray and TeV photons before periastron are mainly produced by the shocked 
electrons around the shock apex and the light curves after periastron are contributed 
by the emission from the shock apex and the shock tail together, which result in the 
asymmetric two-peak light curves. 

Subject headings: binaries: close — gamma rays: stars — pulsars: individual(PSR B 1259-63) 
— X-rays: binaries 
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1. Introduction 

The discovery of the PSR B1259-63/LS 2883 system was first reported in 1992 (Johnston et 
al. 1992), and it is a binary system containing a rapidly rotating pulsar, PSR B 1259-63, in orbit 
around a massive Be star companion LS 2883. The spin period of the pulsar is P = 41.16 ms and 
the spin-down luminosity is L sd - 8 x 10 35 ergs s" 1 . The distance between the system and the 
Earth has been updated to be 2.3 + 0.4 kpc by Negueruela et al. (2011) recently. The emission 
from this system has been widely detected in radio (Johnston et al. 2005), X-rays (Chernyakova 
et al. 2006, 2009; Uchiyama et al. 2009) and TeV gamma-rays (Aharonian et al. 2005, 2009), and 
the light curves are modulated on the orbital period. Especially, the X-ray and TeV light curves are 
similar and display two-peak profiles. Recently, this system was first detected in the GeV range 
by the Fermi satellite in its last periastron passage in 2010 mid-December (Abdo et al. 201 1; Tarn 
et al. 201 1). An interesting GeV flare was observed with a cut-off energy at several hundred MeV, 
which is difficult to explain with the traditional lepton model. 

In the traditional lepton model of gamma-ray binaries hosting a pulsar, the interaction 
between the pulsar wind and the stellar outflow will terminate the winds with a shock roughly at 
the position where the dynamical pressures of the pulsar wind and the stellar wind are in balance, 
and this shock can accelerate electrons to relativistic energies. These accelerated electrons around 
the shock apex will emit broadband nonthermal emission via a synchrotron process for the X-rays 
or an external inverse Compton (EIC) scattering of the thermal photons from the Be star for the 
TeV gamma-rays (Tavani & Arons 1997; Dubus 2006; Khangulyan et al. 2007; Takata & Taam 
2009; Kong et al. 2011). The X-ray light curve reaches a maximum in flux at periastron in 
the simplest models, which is inconsistent with observations. Some authors used some revised 
leptonic models to explain the drop of photon flux towards periastron, for example by introducing 
some non-radiative losses of electrons (Kangulyan et al. 2007) or varying the microphysical 
parameters (Takata & Taam 2009; Kong et al. 201 1). On the other hand, the synchrotron spectrum 
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has a maximum energy no more than 2.36 x 10 8 eV by equating the synchrotron cooling timescale 
with the particle acceleration timescale (see Sect. 2.4), which seems consistent with the Fermi 
observed GeV cut-off energy during the flare. But we should notice that if the acceleration 
efficiency is constant, the maximum energy from synchrotron radiation will not vary along with 
the orbital phase. If this energy corresponds to the cut-off energy in the flaring period, the spectra 
in other periods cannot be explained properly. The EIC process mainly contributes to the photons 
at above 1 GeV, and also cannot explain the observed flare. 

Some previous studies suggested that the GeV flare could be produced by Doppler-boosting 
the synchrotron radiation (Kong et al. 201 1; Tarn et al. 201 1). The interaction between winds in a 
binary system should produce a bow-like structure. Bogovalov et al. (2008, 2012) presented their 
hydrodynamic simulations of the interaction between the relativistic and nonrelativistic winds 
in the PSR B1259-63/LS 2883 system, using both the unmagnetized and magnetized, both the 
isotropic and anisotropic pulsar winds. They found that the bulk motion of the downstream pulsar 
wind electrons can be accelerated from a Lorentz factor ~ 1 around the shock apex to a very large 
Lorentz factor in the bow shock tail. Some previous works (Khangulyan et al. 2008; Dubus, 
Cerutti & Henri 2010) have used the relativistic Doppler-boosting effect to explain the emission in 
gamma-ray binaries. A similar effect should exist in the PSR B1259-63/LS 2883 system. With a 
large bulk Lorentz factor, the emission from the shock tail should be strongly beamed. When the 
line-of-sight is near the beaming direction, we can receive the boosted GeV flux, otherwise, the 
GeV photons disappear. Coincidentally, the true anomaly of the GeV flare (1 10° - 130°) is almost 
the same as the true anomaly corresponding to the direction of the Earth (130°), where the effect 
of Doppler-boosting is the most significant. Tarn et al. (2011) found that the flux of the Fermi 
observations in the flaring period is enhanced by a factor of 5-10, which suggests a Doppler factor 
of around 1.5-2. It is interesting to note that for the X-ray and TeV bands, the second peaks in 
the light curves are also around the GeV flaring period. Therefore these second peaks may be 
produced by the emission from the shock apex and the Doppler-boosted emission from the shock 



-5- 



tail together. 

The two-peak profiles in both the X-ray and TeV light curves are also distinctive features 
of the PSR B1259-63/LS 2883 system. In addition to the extra-contribution from the Doppler- 
boosting effect, some other anisotropic structures in this system may play significant roles on 
this problem. Bogovalov & Khangoulian (2002) has suggested an anisotropic distribution of 
energy flux in the pulsar wind to interpret the torus and jet-like structures in the center of the 
Crab Nebula. We can imagine that the anisotropy of wind is a common phenomenon in pulsars, 
including PSR B 1259-63. If the spin axis of the pulsar is not perpendicular to the orbital plane, as 
the pulsar moving around its companion star, the energy flux injecting to the termination shock 
will be modulated with respect to the orbital phase. This modulation has a two-peak profile and 
further produces asymmetric two-peak profiles in the observed light curves. 

In this paper, we will use the Doppler-boo sting effect to reproduce the flare in GeV light curve 
in detail. We will also use an anisotropic pulsar wind model, together with the Doppler-boo sting 
effect, to explain the asymmetric two-peak profiles in both the X-ray and TeV light curves. 
A variation of the magnetization parameter <x with the distance to the pulsar suggested in our 
previous paper (Kong et al. 201 1) is included in our calculations. The outline of our paper is 
as follows: in Section 2, we introduce our model in detail. We then present our results and the 
comparison with observations in Section 3. Our discussion and conclusion are presented in 
Section 4. 

2. Model Description 

In our model, the broadband emission of the PSR B1259-63/LS 2883 system is mainly from 
the shock-accelerated electrons, both around the shock apex and in the shock tail. Due to the 
interaction between the pulsar wind and the stellar wind, strong shocks will be formed, and the 
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electrons (and positrons) can be accelerated at the shock front of the pulsar wind. This shock will 
also compress the magnetic field in the pulsar wind. The shocked relativistic electrons move in 
the magnetic field and the photon field of the Be star, and emit synchrotron and IC radiation to 
produce the multiband emission. 



2.1. Shock Geometry 

As illustrated in Fig.l, the interaction between the pulsar wind and the stellar wind will form 
a shock with a hollow cone-like structure. The distance from the shock contact discontinuity in 
the shock apex to the pulsar can be determined by 

1/2 

where d is the separation between the pulsar and its companion and r\ is the ratio of the momentum 
fluxes from the pulsar and the massive star. When the pulsar wind is isotropic, the value of rj 
should be L sd /cMv w , where c is the speed of light, M is the mass-loss rate of the massive star and 
v w is the velocity of the stellar wind. The detail of the wind is presented in Kong et al. (2011). 
Away from the apex, the shock surface becomes a hollow cone. The half-opening angle of the 
shock contact discontinuity should be (Eichler & Usov 1993) 

=2/5 

e = 2.i(i --^-W 3 , (2) 

where fj = mm(r], r]~ l ). In this work, because we use an anisotropic pulsar wind model (see Sec. 
2.2), the value of rj in the shock apex should vary in different orbital phases. However, as shown 
by Bogovalov et al. (2012), the effect of the anisotropic pulsar wind on the bow-shock structure is 
relatively moderate, i.e it can not obviously affect the geometry of the bow shock, so we use the 
mean value of < rj > to determine the location of the shock apex and the shape of the bow shock 
in our calculations. The anisotropy in the pulsar wind mainly affects the energy flux injection to 
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the termination shock, i.e. affects the bulk Lorentz factor of the unshocked pulsar wind (see Eq. 
(3)) and the downstream magnetic field (see Eq. (7)), in different orbital phases. 

In our model, we approximate that the observed emission is mainly produced in two regions 
(as illustrated in Fig.l): (1) Region I around the shock apex. The bulk motion of the particle flow 
in this region is assumed non-relativistic, so the radiation is isotropic. (2) Region II in the shock 
tail. As the particle flow propagating away from the shock apex, the bulk Lorentz factor of the 
flow will be increased gradually from r min ^ 1 in this region (Bogovalov et al. 2008, 2012). The 
bulk motion is mildly-relativistic and the emission should be beamed. When the line-of-sight is 
near the beaming direction, as illustrated in Fig. 1, we can receive the Doppler-boosted photons. 
Note that Eq.(2) introduced by Eichler & Usov (1993) is defined at a distance very far from the 
shock apex, where the moving directions of the unshocked pulsar and stellar winds are nearly 
parallel. As shown in Fig. 1, because the moving directions of the unshocked pulsar and stellar 
winds are not the same in Region II, the mean half-opening angle of the particle flow (p should 
be larger than the angle 6 estimated in Eq.(2) and the angle between the direction of the stellar 
photons and the beaming direction should be ~ <p - 9. In this case, even if the line-of-sight is in 
the beaming direction, the directions of photons before and after scattering are not the same, and 
the EIC process will not be completely suppressed. When the particle flow is very far from the 
shock apex and the massive star, it will deviate from the flow direction of Region II due to the 
effects of the orbital motion and the Coriolis forces (Bosch-Ramon & Barkov 2011). So in our 
work, we assume that the electrons beyond Region II do not contribute to the observed flux. The 
bulk Lorentz factor of the flow at the end of Region II is r max . We also assume that the electron 
numbers with different bulk Lorentz factors (from r min 1 to r max ) are the same in calculating the 
synchrotron and EIC radiations. 
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2.2. Anisotropic Pulsar Wind 

As suggested by Bogovalov & Khangoulian (2002), the distribution of energy flux in the 
pulsar wind should be anisotropic and the particle flux can be considered to be more or less 
isotropic. Defining the angle between the pulsar spin axis and the direction from the pulsar to the 
massive star as # PB , the expression of the bulk Lorentz factor of the pulsar wind in the upstream of 
the termination shock as a function of 6 PB should be as following: 

7i = To + 7mSin 2 # PB , (3) 
where y ~ 200, y m ~ 10 6 - 10 7 . Then the isotropic particle flux could be described by 

* = 2 / Sd 2 v (4) 
m e c 2 (y + 3r m ) 

where m e is the rest mass of the electron. As the pulsar orbiting around the massive star, # PB will 
be modulated on the orbital phase (See Fig. 2), and the characteristics of the pulsar wind in the 
termination shock will vary accordingly. 

As shown in Khangulyan et al. (201 la, 201 lb), the EIC cooling will reduce the Lorentz 
factor of the unshocked electrons. We will not consider this effect in our modelling, and the initial 
bulk Lorentz factor parameter of the unshocked pulsar wind y m should be higher if this effect is 
added. Note that because the pulsar wind bow shock is extended, we only obtain an upper limit of 
the wind anisotropy effect. But the two emission regions in our model are only a part of the whole 
bow shock, the effect of the extending bow shock is relatively moderate. On the other hand, we 
can vary the parameter 6 PB to increase the effect of the anisotropic pulsar wind. 

2.3. Magnetization Parameter 

The magnetization parameter <x is defined as the ratio of the magnetic energy density and the 
particle kinetic energy density in the pulsar wind. As estimated in Kong et al. (201 1), by using the 
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pulsar parameters of PSR 1259-63, the magnetization parameter at the light cylinder should be 

^ = 2^W^ ~ 468 X 10 ( 2.5 x W ( 2.3 x 10«cm ) ( 10? ) ( 5.26 x io^ } ' (5) 
where B L is the magnetic field at the light cylinder, r L is the radius of the light cylinder, 
N e± = N m N G] , N m is the e ± multiplicity and N G] ~ 5.26 x 10 31 (5/3 x 10 11 G)(/ , /47.762ms)- 2 s^ 1 
is the Goldreich- Julian particle flow at the light cylinder. In outer gap models (e.g. Cheng, Ho & 
Ruderman 1986a, 1986b; Zhang & Cheng 1997; Takata, Wang & Cheng 2010), the multiplicity 
due to various pair-creation processes could reach 10 4 - 10 5 . Eq. (5) shows that the pulsar wind 
is Poynting-dominated at the light cylinder. In some studies of the Crab Nebula (Kennel & 
Coroniti 1984a, 1984b), the pulsar wind should be kinetic-dominated (cr ~ 0.003) at a distance 
of r s ~ 3 x 10 17 cm from the pulsar. We can imagine that between the light cylinder and the 
termination shock, the magnetic energy will be gradually converted into the particle kinetic 
energy. In our previous work (Kong et al. 2011), we have already shown that a variation of 
magnetization parameter cr with the distance from the pulsar could help us to reproduce the 
two-peak profiles in light curves, and we will use the same variation in this work and describe it as 

Cr = Cr L (-r , (6) 

where r is the distance from the pulsar, and the typical value of the index a is of the order of unity. 

The downstream magnetic field in Region I could be described as (Kennel & Coroniti 1984a, 
1984b) 

«■ J rjjf ■ 7 ' 2 ^4 CD 



2 8cr 2 + 10cr+l [64cr 2 (cr+ l) 2 + 20cr(cr+ 1)+ 1]2 

M 7 = 1 . (8) 

2 16(CT+ 1) 16(CT+ 1) V ' 

Because Region II is far from the shock apex and it is an oblique shock there, the magnetic field 
should be lower than that in Region I. We assume the ratio between the magnetic fields in the two 
emission regions does not change in different orbital phases, and the exact ratio is determined by 
fitting in our work. 
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2.4. Radiation Process 

It is usually assumed that the unshocked cold electron pairs in the pulsar wind can be 
accelerated to a power-law distribution <2(y e ) ~ (y e _ ^T p (7e,min < 7e < 7e,max) in the 
termination shock front and be injected into the downstream post-shock flow, where p is 
the electron distribution index. The minimum Lorentz factor can be determined from the 
conservations of the total electron number L sd /yim e c 2 = J <2(y e )dy e and the total electron energy 
£ S d = f Q(7 e )7 e m e c 2 dy e , and we can acquire y e , min = y x (p - 2) Hp - 1) for p > 2 (Kirk, Ball & 
Skjaeraasen 1999). The maximum Lorentz factor y e ,max 

can be determined by equating the cooling 
timescale of electrons with the particle acceleration timescale as the following form, 



y e , max = J^-i.iixioY' 2 *- 1 ' 2 , (9) 

where e is the electron charge, cr T is the Thompson scattering cross section, £ is the acceleration 
efficiency which is usually less than unity. Hereafter the convention Q x = <2/10 x is adopted for 
the cgs units. 

The electrons will lose their energies through the radiative or adiabatic cooling processes, 
and the evolved electron spectrum n(y e , t) can be obtained from the continuity equation of the 
electron distribution (Ginzburg & Syrovatshii 1964), 

dn(y e ,t) dy e n(y e ,t) . 

— — + — a = 2(7e), (10) 

at dy e 

where y e is the total energy loss rate of the electrons and Q(y e ) is the injection rate. The coefficient 
of the injection rate in Region I rji = Q(y e )/(y e - l)~ p (7 e , m in < 7e < 7e,max) can be calculated 
from Qua = f Q(y e )dy e , where Q tot = L sd /[Am e c 2 (y Q + 2y m /3)] by assuming the typical scale of 
the shock apex is r s (Dubus 2006). We do not know the exact structure and physical conditions in 
Region II, so we assume the electron injection rates are the same in the two regions for simplicity. 
In this case, about half of the pulsar wind electrons are injected to the emission regions. 

Because the cooling and dynamic flow timescales are much smaller than the orbital period in 
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the PSR B1259-63/LS 2883 system, we use dn(y e , t)/dt = to calculate the electron distribution 
at a steady state and acquire (Khangulyan et al. 2007; Zabalza, Parades & Bosch-Ramon 201 1) 



n(y e ) = — 

Ird J re 



Q(y' e W e - (ID 



For a source with the dynamical timescale T dyn , the electron number at a given Lorentz factor that 
can accumulate in the source could be simply calculated by <2(ye)rnin[r c (y e X T dyn] (Moderski et 
al. 2005), where r c is the cooling timescale. Usually it is assumed the dynamical timescale in 
Region I is r dyn = 3r s /c. We also assume the dynamical timescale in Region II is the same with 
that in Region I for simplicity. For the radiative cooling timescale r c , we use the method proposed 
by Moderski et al. (2005) to calculate 

3 ffl Q TJ 

Tc(7e) = -j '—/[l + 7r F KN(7e)], (12) 

4cr T 7 e [/ B U B 

where U B = B 2 /Sn is the magnetic energy density, U» = L sta JAncR 2 is the seed photon energy 
density in Region I where L star is the luminosity of the massive star and R is the distance between 
the emission region and the massive star, F K n - (1 + by 1 - 5 and b = 4y e (2.8kT sS /m e c 2 ) in Region 
I, where k is the Boltzmann constant and T e ff is the effective temperature of the star. In Region II, 
U, and b should be reduced by a factor of D\ and D* respectively, where D* = 1/T(1 -jScos 6*) is 



the Doppler factor, Y is the bulk Lorentz factor of the particle flow, J3 = ^(r 2 - 1)/T and 0* is the 
angle between the direction of the stellar photons and the moving direction of the flow. Note that 
a mono-energetic photon distribution with energy 2.8kT e ff is used here as a good approximation 
of the thermal distribution of the stellar photons (Moderski et al. 2005). 

In our model, the multiband photons from the PSR B1259-63/LS 2883 system are produced 
by the synchrotron radiation and the EIC process of the shock-accelerated electrons. An 
anisotropic inverse-Compton radiation formula is used in our calculations, in which the radiation 
power at frequency v from a single electron with Lorentz factor y e in the comoving frame is given 
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by (Aharonian & Atoyan 1981) 

dv s^krh^,b e ), (13) 

Jy.r 



dg^cosgsc) = 3ot f" v^ TAR 



^ Bl + 2R-fill + R' (14) 
where /? is the Planck constant, £ = hv/(y e m e c 2 ), b e = 2(1 - cos#sc)7e^v s /(m e c 2 ), /zv s <sz hv < 

y e m e c 2 b e /(l + b e ), Bsc is the angle between the injecting photons and the scattered photons, and 

is varied along with the orbital phase. In Region II, this angle in the comoving frame # S c can be 

related to that in the observer frame (p sc by 1 - cos 9 SC = D obs D*(l - cos <^ S c) (Rybicki & Lightman 

1979; Dubus, Gerutti & Henri 2010; Zdziarski et al. 2012), where D obs = l/r(l -/3cos6> obs ) is 

the Doppler factor and # obs is the angle between the line-of- sight and the moving direction of the 

flow, which is modulated on the orbital period. The flux density of the massive star photons / V ^ TAR 

should be n{RJ Riflhv 3 /c 2 [exp(hv/kT eS ) - 1] and n(RJRn) 2 2D 2 hv 3 /c 2 [exp(D*hv/kT eS ) - 1] 

(Rybicki & Lightman 1979; Zdziarski et al. 2012) in Region I and Region II respectively, where 

R t is the radius of massive star, Rj and Ru are the distances between the emission regions and the 

massive star in Region I and Region II respectively. Because we do not know the exact structure 

of Region II, we use a mean value of < Ru > in our calculations. Note that here the massive star is 

assumed to be a black body emitter for simplicity. 

The radiation from the electrons in both the shock apex and the shock tail are included in our 
calculations. The emission in the shock apex is assumed isotropic, and the radiation in the shock 
tail is beamed. For a point-like source, the Doppler-boo sting effect will increase the photon energy 
by a factor of D obs and increase the detected flux by a factor of D 3 oh& (Dubus, Gerutti & Henri 
2010). We treat Region II as a ring-like shape on the cone-like termination shock approximation 
and we integrate over the different parts of Region II to get a more accurate result. According to 
Eq.(9), the synchrotron spectrum will cut off at the energy of 

3hey 2 B 

/zVsyn,max(#obs) = A)bs(flobs) A ~ 2.36 X lO^D obs (0 obs )eV. (15) 
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The variation of D obs will help us to obtain the observed cut-off energy in the flaring period and 
the flare structure in the light curve. 

3. Results 

In this section, we will present some calculated results using our model, and compare them 
with the observations. We use the updated parameters for the PSR B1259-63/LS 2883 system 
in our calculations as follows (Negueruela et al. 2011): For the orbital parameters, we take the 
eccentricity e = 0.87, the semimajor axis a = 7 AU; For the compact object PSR B 1259-63, 
we take the spin-down luminosity L sd = 8 x 10 35 erg s _1 ; For the Be star LS 2883, we take the 
stellar luminosity L star = 7.3 x 1O 4 L , the stellar radius R t = 10R Q , the effective temperature of 
the star T eS = 30000 K. The distance between the system and the Earth is taken to be 2.3 kpc. 
The angle between the line-of-sight and the orbital plane is taken to be 65°, and the true anomaly 
corresponding to the direction of the Earth is taken to be 130°. The input parameters we used 
are as follows: The mean ratio of the momentum fluxes from the pulsar and the massive star is 
< r] >= 0.16, which corresponds to the mass-loss rate of the massive star M - 2.6 x 10~ 8 M yr _1 
with a velocity of v w ~ 10 8 cm s _1 , which is consistent with the typical value of the polar wind 
in the Be star (Waters et al. 1988). The corresponding half-opening angle of the shock cone is 
6 ~ 58°. We choose the half-opening angle of Region II <p = 65°, so that around the true anomaly 
130°, we can receive strong beamed emission from the shock tail. The angle between the direction 
of stellar photons and the moving direction of the flow in Region II is taken as 8 t ~ (p - 6. The 
maximum bulk Lorentz factor contributing to the observed flux in the shock tail is r max = 2.0, the 
bulk Lorentz factor parameters of the unshocked pulsar wind are y = 200 and y m = 2.0 x 10 6 , the 
true anomaly of the projection of the pulsar spin axis in the orbital plane is 20°, the angle between 
the pulsar spin axis and the orbital plane is 46°, the magnetization parameter at the light cylinder 
is cr L = 8 x 10 3 (for N m ~ 5.8 x 10 4 ) and the decay index is a = 1.1, the electron distribution index 
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is /? = 2.1, the acceleration efficiency is £ = 0.36. The ratio of magnetic fields between Region 
II and Region I is assumed to be 0. 1 . The mean distance between the massive star and Region 
II < Ru > is taken to be 2.0 times of the binary separation, which is consistent to the results in 
Bogovalov et al. (2008, 2012) for r max = 2.0. Note that by using the parameters introduced above, 
the dynamical timescale Td yn in Region I is ~ d/c and that in Region II is less than ~ 2d/c. The 
difference between these two timescales is small, so the assumption that these two regions have 
the same dynamical timescales in our calculations is reasonable. The periastron is taken to be at 
orbital phase ~ throughout the paper. 

In Fig. 2, we show the variations of the angle between the pulsar spin axis and the line joining 
the two stars 9 BP , the bulk Lorentz factor of the unshocked pulsar wind ju the magnetization 
parameter at the termination shock cr, and the magnetic field around the shock apex B with respect 
to the orbital phase. We can see that because the pulsar spin axis is not perpendicular to the orbital 
plane, # B p is modulated on the orbital phase. For the case that the true anomaly of the projection 
of the pulsar spin axis in the orbital plane is 20° and the angle between the pulsar spin axis and the 
orbital plane is 46°, #bp will reach the minimum and maximum of 46° and 134° at the true anomaly 
of 20° and -160° respectively. By using the anisotropic energy flux injection of the pulsar wind 
described in Section 2.2, j\ will vary with respect to # B p (See Eq. (3)) and reach the maximum 
when 6 BP = 90° at the true anomaly of -70° and 110°. In this case a two-peak profile appears in 
the distribution of j\ . This two-peak structure will affect the minimum Lorentz factor y e>m i n in 
the electron distribution, and further affects the electron number at a certain Lorentz factor. The 
two-peak distribution on the electron number at a certain Lorentz factor will help us to obtain 
the two-peak profiles in light curves. We can also see that the magnetization parameter at the 
termination shock cr is modulated with respect to orbital phase, which is due to the different shock 
distance r s in different orbital phases. The magnetic field in the shock apex B is also modulated 
with respect to orbital phase, which is similar to the modulation of magnetization parameter cr. 
But this modulation is not symmetrical because of the anisotropy of the pulsar wind (See Eq. (3) 
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and Eq. (7)). 

Here we use the parameters introduced above to calculate the spectra and the comparisons 
with observations are shown in Fig. 3. The related timescales are shown in Fig. 4. It can be seen 
that our calculated X-ray and GeV emission is mainly produced by the synchrotron radiation and 
the TeV photons are mainly contributed by the EIC effect where the relativistic shocked electrons 
up-scatter the soft photons from the massive star. The X-ray and TeV flux at true anomaly of 
120° is contributed by the emission from shock apex and shock tail together, but the flux at 
true anomaly of -60° is mainly produced by the electrons in the shock apex. By choosing the 
acceleration parameter £ = 0.36 in our modelling, the cut-off energy of the synchrotron spectrum 
in the unboosted region is hv syn ^ max ~ 84 MeV, which is a little lower than the lower limit of the 
energy range of the observations (100 MeV). So in the lower panel of Fig. 3, the pre-periastron 
synchrotron spectrum only have a small contribution in the above 100 MeV range. In the upper 
panel of Fig. 3, the line-of-sight is near the beaming direction. The photon flux and the photon 
energy will be strongly Doppler-boosted, and the cut-off energy in the synchrotron spectrum could 
be boosted to /zv synimax ~ 300 MeV, consistent with the Fermi observed GeV cut-off energy during 
the flare (Abdo et al. 201 1; Tarn et al. 201 1). In this case, the observed flare can be reproduced. 
Unlike in the synchrotron component above 100 MeV, the EIC component above 1 TeV is not 
dominated by the boosted emission from the shock tail in the upper panel of Fig. 3. This is 
because (1) in the comoving frame the seed photon density in the shock tail is lower than that 
in the shock apex reduced by the Doppler de-boosting; (2) because of the anisotropy of the EIC 
emission, the flux will be suppressed when the angle between the directions of the input and output 
photons is small. In Fig. 4, we can see that the radiative cooling of the electrons with Lorentz 
factor 10 5 - 10 6 is mainly dominated by the EIC radiation in the Klein-Nishina (KN) regime, 
and the radiative cooling of the electrons with higher energies is dominated by the synchrotron 
radiation. In our calculations, the minimum Lorentz factor of the shocked electrons is ~ 10 5 yi,6, 
so the EIC process will never be in the Thomson regime. As a result, the synchrotron component 
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is higher than the EIC component in our calculated spectra. Our calculated synchrotron spectrum 
in the lower panel of Fig. 3 does not fit the observations in the GeV range well, but note that the 
spectrum data in pre-periastron period given by different groups (Abdo et al. 2011; Tarn et al. 
201 1) also are not consistent with each other, which may be due to the low value of the photon 
flux. 

Our calculated broadband light curves of the PSR B1259-63/LS 2883 system and the 
comparisons with observations are presented in Fig. 5. We can see that the asymmetric two-peak 
profiles in the X-ray and TeV light curves can be well reproduced by our model. The observations 
before periastron are mainly contributed by the emission in the shock apex. The shock apex region 
also contributes to some emission in the post-periastron range, and the whole light curves display 
the asymmetric two-peak profiles with the flux before periastron higher than that after periastron. 
The anisotropic pulsar wind plays a significant role in producing this asymmetric two-peak profile 
because of its anisotropic energy flux injection to the termination shock. The emission from the 
shock tail is unimportant before periastron, because the angle between the line-of-sight and the 
direction of the flow in the shock tail O ^ S is large. But it has an obvious contribution in the light 
curves after periastron, i.e. the second peaks in both the X-ray and TeV light curves are produced 
by the electrons in the shock apex and the shock tail together. For the GeV light curves, we can 
see that almost all the flux in the GeV flares are produced by the emission from the shock tail, 
where the emission is beamed and both the photon energy and detected flux are Doppler-boosted. 



4. Conclusions and Discussion 

The PSR B1259-63/LS 2883 system is an attractive binary system and special for having 
distinct modulations in the light curves of different energy bands. The X-ray and TeV light curves 
are similar and have two peaks before and after periastron respectively. The GeV light curve has 
an interesting flare several weeks after the periastron passage, whose spectrum cuts off at several 
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hundred MeV. In this paper, we have modelled the X-ray, GeV and TeV observations of this 
system. In our calculations, an anisotropic pulsar wind model and the relativistic Doppler-boosting 
effect are considered, and the X-ray and GeV photons are mainly from the synchrotron radiation 
and the TeV emission is mainly from the EIC process. We also assume that the emission from 
both the shock apex and the shock tail can contribute to the observations, and the shock apex 
radiation is isotropic and the emission from the shock tail is beamed because of the flow in the tail 
moving mildly-relativistically. We find that for the X-rays and TeV gamma-rays, the anisotropic 
energy flux injection of the pulsar wind plays an important role in producing the asymmetric 
two-peak profiles in the light curves. The observations before the periastron are mainly produced 
by the shocked electrons around the shock apex and the photons after periastron are contributed 
by the emission from the shock apex and the shock tail together. For the GeV gamma-rays, 
the observed flare is contributed by the emission from the shock tail, where the synchrotron 
photons are Doppler-boosted strongly. Unfortunately, we do not know the exact structure and 
physical processes in the shock tail, so in our calculations we make some assumptions, and we 
just show that our model is a possible way to reproduce the multiwavelength features in the 
PSR B1259-63/LS 2883 system. A more exact modelling could be done by using the detailed 
simulation results on binary pulsar systems (Bogovalov et al. 2008, 2012; Takata et al. 2012). 

In our modelling, we use an assumption that the magnetization parameter a varies with the 
distance to the pulsar, which is in principle possible. Some recent studies found that the variation 
of magnetization parameter cr could be in a different way. For example, Aharonian, Bogovalov 
& Khangulyan (2012) showed that the magnetization parameter <x should decrease abruptly 
within 10 10 cm from the pulsar by fitting the observations of the Crab pulsar. Because the exact 
conversion process from Poynting flux to kinetic energy in the pulsar wind is still unclear, we think 
both models cannot be excluded. Some previous models implied that the energy conversion could 
exist over the entire distance from the pulsar to the termination shock (Coroniti 1990; Lyubarsky 
& Kirk 2001), and Contopoulos & Kazanas (2002) further showed that <x could decrease inversely 



-18- 



proportional to the distance from the light cylinder. In our previous paper (Kong et al. 2011), a 
variation of the magnetization parameter cr with the distance from the pulsar could help us to 
reproduce the two-peak profiles in light curves, and we also need a variation index a = 1 . 1 to fit 
the observations better in this work. The adoption of a = 1.1 here is purely phenomenological, 
which is affected by both the unclear magnetic energy dissipation and the EIC cooling of the 
unshocked pulsar wind (Khangulyan et al. 201 la, 201 lb). 

In our previous paper (Kong et al. 2011), we discussed the effect of the disk in the stellar 
wind on the X-ray and TeV light curves. The existence of a disk in the PSR B1259-63/LS 2883 
system was confirmed by the radio observations (Johnston et al. 1996, 2005), but the exact 
position of the disk is still unclear. Some radio observations suggested that the disk is tilted with 
respect to the orbital plane and the line of intersection between the disk plane and the orbital plane 
is oriented at about 90° with respect to the major axis of the binary orbit (Wex et al. 1 998; Wang, 
Johnston & Manchester 2004). Chernyakova et al. (2006) further suggested that the half-opening 
angle of the disk (projected on the pulsar orbital plane) is A# disk 18°. 5, and the intersection 
between the stellar equatorial plane and the orbital plane is inclined at #disk - 70° to the major 
axis of the pulsar orbit by fitting the X-ray and TeV light curves. Our previous calculations (Kong 
et al. 201 1) showed that the X-ray flux increases in the passage of the disk, but the flux in the 
TeV range decreases significantly, which is consistent with the analysis by Kerschhaggl (2011). 
But note that we only consider the emission from the shock apex in Kong et al. (2011). When 
the pulsar entering the disk, the mass flux density will increase by a factor of 30-100 (Waters et 
al. 1988). By assuming the stellar wind velocity reduced by a factor of 10, the momentum flux 
density ratio rj will decrease by a factor of 3-10 in the disk. In this case, the half-opening angle of 
the emission region in the shock tail will decrease and the angle between the line-of-sight and the 
direction of the flow in the shock tail 6 obs will increase. The effect of Doppler-boosting will be 
suppressed and the emission from shock tail will be reduced. We can see from Fig. 5 that there 
are only upper limits in GeV light curves between true anomaly 90° and 115°. The non-detection 
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in this period may be due to the suppression of the shock tail emission in the disk. The TeV 
observations also show a reduction in flux during the estimated disk passage, which are consistent 
with the explanation that emission from both the shock apex (Kong et al. 201 1) and the shock tail 
will decrease in the disk. In X-ray band, there are significant signals between true anomaly 90° 
and 1 15°. But note that the X-ray flux from the shock apex could increase in the passage of the 
disk (Kong et al. 2011). The increase of the shock apex emission and the decrease of the shock 
tail emission in the disk will compete with each other, and the total flux could have no significant 
change. 

In a recent paper, Khangulyan et al. (2011a) investigated the emission spectrum produced 
by the EIC process of the unshocked pulsar wind. They argued that the gamma-ray flare in GeV 
band after the periastron can be explained by the EIC emission of the cold pulsar wind with the 
bulk Lorentz factor y « 10 4 . The Be star disk plays a significant role in producing GeV flare in 
their model. First, the radiation of the shocked stellar disk can provide a dense photon target for 
the EIC scattering (van Soelen & Meintjes 201 1). Second, the strong ram pressure inside the disk 
makes the wind termination shock stand close to the pulsar, and then the EIC luminosity should 
be suppressed. When the pulsar escapes the disk, the unshocked pulsar wind zone towards the 
observer is significantly increased. Consequently, an enhancement of the gamma-rays will be 
observed. Although in our model the GeV flare is produced by the Doppler-boosted synchrotron 
radiation, the EIC process of the unshocked pulsar wind may also have some contributions. 
Note that in the pre-periastron spectrum data in Tarn et al. (2011), the Fermi detected emission 
is concentrated in a narrow band between 1-25 GeV. The observed spectrum is different from 
the spectrum produced by synchrotron radiation, and is similar to the calculated spectra in 
Khangulyan et al. (201 la). If the data in this range is indeed correct, we argue that this part of 
GeV photons may be produced by the EIC process of the unshocked electrons before termination 
shock in the pulsar wind. 
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If the 3D structure of the shock is a simple hollow cone at the shock tail, we expect that 
two flares per orbital period will be observed when the angle between the line-of-sight and the 
orbital plane is smaller than the half-opening angle of the shock cone, and each flare can last for 
about 2/r orbital phase ~ 12°(ryi0) -1 , where T is the bulk Lorentz factor at the tail of the shock. 
Otherwise, when the line-of-sight is near the edge of the hollow cone or outside the hollow cone, 
as illustrated in Fig. 1, only one flare could be observed in the light curve. As the angle between 
the line-of-sight and the flow direction increasing, the flare will be smoother and disappear 
eventually. In principle, our model could be used in other similar gamma-ray binary systems in 
our Galaxy with different observational angles and distinct light curves. 
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Fig. 1. — Geometry of the termination shock. The interaction between the pulsar wind and the 
stellar wind forms a termination shock with a hollow cone-like structure. Region I is around the 
shock apex, and the particle flow there is moving non-relativistically and radiating isotropically. 
Region II is in the shock tail, and the particle flow there is moving mildly-relativistically and the 
emission is beamed. 
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Fig. 2. — From top to bottom: the variations of the angle between the pulsar spin axis and the line 
joining the two stars # B p, the bulk Lorentz factor of the unshocked pulsar wind y\, the magneti- 
zation parameter at the termination shock cr and the magnetic field around the shock apex B with 
respect to the orbital phase. 
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Fig. 3. — The calculated spectra as compared with observations. The spectrum data are taken 
from Aharonian et al. (2005), Abdo et al. (2011) and Tarn et al. (2011). In the upper panel, the 
spectra are calculated at true anomaly of 120° and the data are taken in the post-periastron range. 
In the lower panel, the spectra are calculated at true anomaly of -60° and the data are taken in the 
pre-periastron range. The dashed lines and dotted lines correspond to the emission from the shock 
apex and shock tail respectively, and the solid lines correspond to the total flux. 
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Fig. 4. — The calculated timescales with respect to the electron energies. The solid lines corre- 
spond to the dynamical timescale r dyn , the dashed lines and dotted lines correspond to the cooling 
timescales r c in the shock apex and shock tail respectively. In the upper panel, the timescales are 
calculated at true anomaly of 120°, and in the lower panel, the timescales are calculated at true 
anomaly of -60°. 
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Fig. 5. — The calculated multi-band light curves as compared with observations. The X-ray light 
curve data are taken from Chernyakova et al. (2006, 2009; circles) and Abdo et al. (2011; trian- 
gles), the > 100 MeV light curve data are taken from Abdo et al. (201 1), the > 200 MeV light curve 
data are taken from Tarn et al. (201 1) and the TeV light curve data are taken from Aharonian et al. 
(2005, 2009). The dashed lines and dotted lines correspond to the emission from the shock apex 
and shock tail respectively, and the solid lines correspond to the total flux. The empty triangles are 
upper limits in observations. The vertical dashed lines correspond to the estimated disk passage. 



